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Abstract 

In  this  paper  we  present  a  hierarchy  of  models  built  to  describe  the  overall  performance  of  a  single  H2S  fuelled  button  cell  solid  oxide  fuel 
cell  (SOFC).  The  cell,  used  in  the  experimental  studies  of  Liu  et  al.  [M.  Liu,  G.  Wei,  J.  Luo,  A.R.  Sanger,  K.T.  Chuang,  Use  of  metal  sulfides  as 
anode  catalysts  in  H2S-air  SOFCs,  J.  Electrochem.  Soc.  150  (2003)  1025-1029],  was  a  planar  cell  with  a  circular  disc-like  electrode  assembly 
and  the  fuel  and  air  flowing  through  a  concentric  cylindrical  tube  assembly.  The  goal  is  to  model  the  electrochemical  reaction  coupled  with  mass 
transfer,  fluid  flow  and  current/voltage  distribution  in  an  yttria  stabilized  zirconia  electrolyte  fuel  cell  assembly  operated  between  750  and  850  °C. 
The  models  built  range  in  complexity  from  an  algebraic  system  of  equations  that  calculates  the  activation,  concentration  and  ohmic  losses,  to  a 
two-dimensional  finite  element  model  that  solves  all  the  physics  in  the  SOFC  simultaneously.  Kinetic  parameters  in  these  (progressively  more 
comprehensive)  models  have  been  estimated  and  compared,  leading  hopefully  to  more  accurate  estimates  for  these  parameters. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Hydrogen  sulfide  (H2S)  is  a  by-product  of  the  natural  gas 
and  the  petrochemical  industries.  In  addition  to  being  toxic  to 
most  forms  of  life  it  is  an  unpleasant  smelling  air  pollutant.  The 
primary  industrial  process  used  to  dispose  off  H2S  is  the  Claus 
process  [1].  The  Claus  process  converts  the  gaseous  H2S  to  sul¬ 
fur,  a  solid  at  ordinary  temperatures,  using  a  fairly  complicated 
multistage  process  that  has  high  operating  costs. 

If  H2S  can  be  used  as  fuel  in  a  fuel  cell,  it  can  be  disposed  off 
in  a  much  simpler  process  that  generates  high  quality  electrical 
energy.  The  earliest  studies  on  the  possibility  of  using  H2S  in  a 
fuel  cell  were  conducted  by  Pujare  et  al.  [2,3] .  Since  then  a  num¬ 
ber  of  studies  have  been  published  on  fuel  cells  running  on  H2S 
[4-12].  These  studies  were  experimental  ones  that  evaluated  dif¬ 
ferent  materials  for  the  anode,  anode  catalysts,  and  electrolyte. 
While  some  work  has  been  reported  on  using  H2S  in  fuel  cells 
with  proton  conducting  electrolytes  [9-12],  most  studies  have 
used  oxide  ion  conducting  electrolytes. 

For  SOFCs  using  H2S  as  fuel,  many  different  anode  materi¬ 
als  have  been  investigated  including  thiospinels  [2,3],  Pt  [13,5], 
various  metal  sulfides  [4,6,8],  and  lanthanum  strontium  vana¬ 


*  Corresponding  author.  Tel:  +1  780  492  5810;  fax:  +1  780  492  2881. 
E-mail  address:  kumar.nandakumar@ualberta.ca  (K.  Nandakumar). 

0378-7753/$  -  see  front  matter  ©  2006  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2006.07.014 


date  (LSV)  [7].  Although  all  of  these  studies  present  fuel  cell 
performance  data  in  the  form  of  current  density  and  or  power 
density  curves,  the  main  goal  of  these  studies  was  to  find  mate¬ 
rials  that  lead  to  stable  performance  for  H2S  fuel  cells  without 
degradation  over  time. 

Among  the  more  recent  studies,  a  group  at  the  University  of 
Alberta,  Liu  et  al.  [6]  and  Wei  et  al.  [8,14]  have  presented  per¬ 
formance  data  for  H2S  fuelled  SOFCs  using  mixed  sulfide  (and 
mixed  sulfide,  metal  and  electrolyte  composite)  anodes.  These 
data  include  current  density  curves  at  different  temperatures  and 
fuel  and  air  flow-rates.  Although  some  impedance  data  are  given, 
these  data  are  for  complete  cells  and  does  not  attempt  to  iso¬ 
late  anode  processes  from  cathode  processes.  Another  recent 
study  from  the  Georgia  Institute  of  Technology  by  Aguilar  et  al. 
[7]  presents  performance  data  at  different  temperatures  from  a 
SOFC  using  an  LSV  anode.  This  work  also  gives  performance 
and  impedance  data  at  different  fuel  compositions  to  demon¬ 
strate  the  selectivity  of  the  anode  to  H2S  in  preference  to  H2. 
Neither  of  these  studies,  however,  pin  down  the  reaction  mech¬ 
anism  or  give  any  kinetic  parameters. 

It  is  in  the  presence  of  such  limited  data  and  uncertainties 
about  detailed  mechanisms  that  we  attempt  to  develop  a  math¬ 
ematical  modelling  framework.  Such  models  can  still  be  useful 
in  assessing  the  relative  importance  of  transport  versus  reaction 
mechanisms. 
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Nomenclature 

ak 

activity  of  species  k 

61k,  in 

activity  of  species  k  at  inlet 

ci 

integration  constant  for  voltage  profile  in  phase  j 

ck 

molar  concentration  of  k  (molem-3) 

Ct 

total  molar  concentration  (molem-3) 

Deff,k 

effective  diffusivity  of  component  k  (m2  s-1) 

Dij 

multicomponent  composition  dependent  diffusiv¬ 
ity  of  species  i  in  species  j  (m2  s-1) 

Dik 

multicomponent  Maxwell-Stefan  diffusivity  of 
species  i  in  k  (m2  s-1) 

Dfc,k 

Knudsen  diffusivity  for  k  (m2  s-1) 

Eoc 

open  circuit  voltage  (V) 

£° 

open  circuit  voltage  for  unit  activities  of  all  par¬ 
ticipating  species  (V) 

E0' 

formal  potential  (V) 

f 

function  representing  the  1-D  explicit  model 

F 

Faraday’s  constant  (96,487  C  mole-1) 

g 

function  representing  the  1-D  implicit  and  2-D 
models 

i 

cell  current  density  (Am-2) 

exchange  current  density  for  electrode  j  (A  m-2) 

'/.1 

limiting  current  density  for  electrode  j  (A  m-2) 

lm 

current  density  calculated  using  the  1-D  implicit 
and  2-D  models  (Am-2) 

^H20 

^H2S,in(2F  -  FH20,in)/FH20,in(2F  +  FH2S,in) 

Ks2 

^H2S,in(^  -  Ps2M)/ Ps2M(?P  +  FH2S,in) 

h 

length  or  thickness  of  j  (m) 

mj 

integration  constant  for  voltage  profile  in  phase  j 

n 

number  of  i-V  data  points  used  in  parameter  esti¬ 
mation  routine 

ne 

number  of  electrons  transferred 

P 

gauge  pressure,  solved  for  in  the  flow  sub- 
domains  in  the  2-D  models  (Pa) 

P 

total  pressure  (bar) 

Pk 

partial  pressure  of  k  (bar) 

Pk,m 

inlet  partial  pressure  of  k  (bar) 

r 

distance  in  radial  direction  (m) 

R 

gas  constant  (8.314  Jmole-1  K-1) 

^contact 

sum  of  contact  and  lead  ohmic  resistances  (£2) 

T 

temperature  (K) 

u 

radial  component  of  mass  average  velocity 
(ms-1) 

V 

axial  component  of  mass  average  velocity  (ms-1) 

V 

molar  average  velocity  (ms-1) 

V 

temperature  corrected  volumetric  flow-rate 
(m3  s-1) 

Fcell 

cell  voltage  (V) 

Vm 

cell  voltage  calculated  using  the  1-D  explicit 
model  (V) 

Oik 

mass  fraction  of  component  k 

Xk 

mole  fraction  of  component  k 

z 

distance  in  axial  direction  (m) 

Vectors  and  matrices 
I  identity  matrix 

i  current  density  vector  (A  m-2) 

ji  diffusive  mass  flux  of  component  i  (kg  m-2  s-1) 

n  normal  vector  (m) 

the  total  mass  flux  of  the  kth  component 
(kgm-2  s-1) 

p,  p*  vectors  of  parameter  values  needed  to  solve  the 
model  equations 
t  tangent  vector  (m) 

v  mass  average  velocity  vector  (ms-1) 

Greek  letters 

ft  charge  transfer  coefficient  in  Butler-Volmer 

equation 

€j  porosity  of  phase  j 

?7act  activation  losses  (V) 

?7conc  concentration  losses  (V) 

kj  permeability  in  electrode  j  (m2) 

/ 1  viscosity  (kg  (ms)-1) 

cpj  electrical  potential  of  phase  j 

p  density  (kgm-3) 

oj  electrical  conductivity  of  phase  j  (Q~l  m-1) 

Tj  tortuousity  of  phase  j 

A stoichiometric  coefficient  for  component  k  equi¬ 
librium  potential  for  electrode  j  (V) 


There  are  two  possible  overall  reactions  that  use  H2S  directly 
in  a  SOFC: 

H2S  +  j02  -vr-  H20  +  jS2  (1) 

H2S  +  |02  O  H20  +  S02  (2) 

At  the  operating  temperatures  of  a  SOFC,  these  reactions  are 
equally  favoured  thermodynamically  and  thus  produce  equilib¬ 
rium  voltages  that  are  very  close  to  each  other,  e.g.,  0.782  V  for 
the  first  and  0.75  V  for  the  second  at  800  °C  [8].  Besides  these 
reactions,  part  of  the  H2S  fed  dissociates  into  H2  and  S2  and  the 
H2  can  then  undergo  electrochemical  oxidation. 

It  is  not  yet  fully  clear  which  (if  any)  of  these  reactions 
predominates  in  a  SOFC  running  on  H2S.  Some  of  the  above 
mentioned  studies  have  attempted  to  figure  out  which  of  the 
above  electrochemical  reactions  predominates  on  the  anode  by 
looking  at  the  open  circuit  voltage  (OCV)  [8]  and  analyzing  the 
anode  exhaust  gases  [4].  However,  the  only  broadly  accepted 
conclusion  drawn  has  been  that  a  H2S  SOFC  produces  more 
SO2  at  higher  loads  and  more  S2  closer  to  OCV  [7]. 

1.1.  Modelling  SOFCs 

Mathematical  models  of  fuel  cells  are  design  and  optimiza¬ 
tion  tools  that  can  help  in  selecting  operating  conditions,  cell 
materials,  stack  design,  controller  design,  among  other  things. 
The  above  capabilities  depend  on  the  ability  to  quickly  mod¬ 
ify  operating  conditions  or  revamp  cell  design  in  the  model  and 
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examine  the  effects  on  cell  or  stack  performance.  Thus,  models 
enable  simulations  or  numerical  experiments  that  can  dramati¬ 
cally  decrease  development  time  for  a  fuel  cell  design  [15]. 

As  with  any  fuel  cell  system,  SOFCs  can  be  modelled  on 
many  levels.  Models  have  been  built  at  the  electrode,  single 
cell,  and  cell  stack  levels.  The  earliest  studies  by  Debendetti 
and  Vayenas  [16]  were  on  the  single  cell  level.  This  model  was 
later  extended  to  model  a  layer  of  cells  in  a  cross-flow  stack  [17]. 
These  models  were  2-D  in  nature  and  only  modelled  the  average 
concentrations  and  temperature  in  each  cell  in  the  stack  layer, 
while  ignoring  mass  transfer  and  electrochemical  kinetics.  Since 
these  studies,  a  large  number  of  models  with  varying  levels  of 
complexity  and  detail  have  been  built,  and  a  good  review  can  be 
found  in  ref.  [18]. 

In  this  work  we  are  modelling  the  fuel  cell  investigated  by 
Liu  et  al.  [6].  The  electrodes  and  electrolyte  assembly  is  in  the 
shape  of  a  circular  disc  held  between  the  air  outlet  tube  on  one 
side  and  the  fuel  outlet  tube  on  the  other.  On  both  sides,  inflow 
occurs  through  the  inner  one  and  the  outflow  through  the  outer 
tube.  Both  the  air  and  fuel  sides  are  at  atmospheric  pressure. 
This  geometry  is  very  common  and  is  also  known  as  a  “button 
cell”.  The  data  sets  that  are  used  in  this  work  to  validate  the 
models  are  from  a  sample  cell  that  used  Co-Mo-S  mixed  with 
Ag  as  anode,  and  Pt  as  cathode.  The  electrolyte  was  8  mole% 
YSZ.  The  cell  assembly  and  the  gas  feeds  and  outlets  are  held 
in  a  temperature  controlled  furnace. 

In  this  paper  we  present  four  different  models  of  increas¬ 
ing  complexity  and  examine  their  fit  to  experimental  data  given 
in  Ref.  [6].  Of  these,  the  first  two  are  one-dimensional  models 
that  solve  for  the  electrochemistry  as  well  as  electrical  potential 
and  mass  transfer  in  the  electrodes.  We  start  with  the  simplest 
model  possible  that  explicitly  calculates  the  output  voltage  of 
the  cell  given  the  current  being  drawn  from  it.  This  model  is 
derived  using  a  one-dimensional  analysis  along  the  axis  of  the 
cell  and  calculates  the  three  voltage  loss  terms:  the  electrical 
resistance  losses,  the  mass  transfer  losses,  and  the  electrochemi¬ 
cal  activation  losses  separately  by  ignoring  the  coupling  between 
the  mass  transfer  and  electrochemistry.  The  second  model  is  a 
refined  version  of  the  first  where  the  mass  transfer  losses  and 
the  electrochemical  activation  losses  are  coupled  and  cannot  be 
calculated  explicitly.  For  this  model,  the  performance  of  the  cell 
is  calculated  by  solving  a  set  of  coupled  non-linear  algebraic 
equations.  The  third  and  fourth  models  are  two-dimensional  axi- 
symmetric  models  where  the  flow,  mass  transfer,  voltage  fields, 
and  electrochemical  reaction  are  all  modelled.  The  fourth  model 
is  identical  to  the  third,  except  the  diffusive  mass  transfer  is  mod¬ 
elled  using  the  multi-component  Maxwell-Stefan  formulation 
instead  of  Fick’s  law.  The  resulting  coupled  partial  differential 
equations  are  solved  using  finite  element  methods  (FEM). 

In  order  to  solve  these  models,  various  parameters  such  as 
diffusion  coefficients,  viscosities,  electrical  conductivities  of 
the  electrodes  and  electrolyte,  and  electrochemical  reaction  rate 
parameters  need  to  be  defined.  While  many  of  these  parameters 
are  known  or  can  be  estimated,  some,  such  as  exchange  current 
densities  for  the  electrochemical  reactions  in  the  electrodes  are 
not  known  with  any  certainty.  There  are  no  data  on  electrochem¬ 
ical  kinetic  parameters  for  H2S  anodes  in  literature.  While  there 


are  some  kinetic  data  for  O2  reduction  on  Pt  cathodes,  there  is  lit¬ 
tle  agreement  in  the  SOFC  research  community  on  the  reaction 
mechanisms  [19].  Since  there  is  no  reliable  data  for  the  elec¬ 
trochemical  kinetic  parameters,  we  use  non-linear  least  squares 
estimation  in  this  work  to  estimate  the  kinetic  parameters  in  all 
four  models. 

In  the  body  of  this  paper,  we  start  with  a  detailed  development 
of  all  the  models  used  and  discuss  the  model  parameters  needed. 
We  then  introduce  the  mathematical  formulation  used  for  esti¬ 
mating  the  unknown  parameters.  The  results  of  the  parameter 
estimation  are  presented  next,  and  we  examine  these  estimates 
and  the  quality  of  the  fit  of  the  models  to  the  experimental  data 
to  draw  conclusions  on  the  fidelity  and  usefulness  of  the  mod¬ 
els.  We  end  with  some  thoughts  on  how  these  models  can  be 
improved  and  plans  for  future  work. 

2.  Modelling 

The  1-D  and  2-D  model  geometries  are  given  in  Fig.  1.  In 
developing  the  2-D  models,  the  cylindrical  geometry  of  the 
experimental  cell  and  the  symmetry  around  the  central  axis  is 
used  to  reduce  the  modelling  domain  to  a  2-D  axi-symmetric 
domain  as  shown  in  the  top  section  of  Fig.  1.  The  sub-domains 
in  this  reduced  domain  are  (going  left  to  right):  (i)  the  fuel  flow 
channel,  £2/?fuei  (ii)  the  anode,  Qa  (iii)  the  electrolyte,  Qm  (iv)  the 
cathode,  Qc  and  (v)  the  air  flow  channel,  %air-  The  1-D  model 
geometry,  given  in  the  bottom  section  of  Fig.  1  does  not  con¬ 
sider  the  flow  channels  as  the  inherently  2-D  flow  field  cannot 
be  modelled  correctly  in  a  1-D  model. 

In  the  2-D  models,  the  current  collectors  are  located  at  the 
vertical  interface  between  the  flow  and  electrode  sub-domains 
on  both  sides.  In  the  1-D  model,  the  anode  current  collector  is  the 
left- mo st  point  of  the  domain  while  the  cathode  current  collector 
is  the  right-most  point  in  the  domain. 

The  main  assumptions  used  in  the  models  presented  are: 

(1)  All  models  presented  are  isothermal. 

(2)  The  only  electrochemical  reaction  modelled  on  the  anode 
is  the  oxidation  of  H2S  to  H2O  and  S2  (reaction  (14))  even 
though  there  are  probably  other  reactions  occuring  simulta¬ 
neously. 

(3)  The  open  circuit  voltage  is  taken  as  the  zero  current  voltage 
in  the  data  set  being  modelled. 

(4)  In  the  1  -D  models,  the  absolute  pressure  is  assumed  constant 
in  the  electrodes. 

We  explain  these  assumptions  in  more  detail  in  relevant  sub¬ 
sections. 

The  assumption  that  all  models  are  isothermal  is  justified  in 
the  current  study  for  the  following  reasons: 

(1)  The  experimental  cell  modelled  is  in  a  temperature  con¬ 
trolled  electric  furnace. 

(2)  The  current  density  drawn  from  the  cell  is  relatively  small, 
which  means  the  heat  generation  is  also  small. 

(3)  Although  the  fuel  and  air  streams  are  not  preheated  before 
they  are  fed  to  the  SOFC,  the  flow-rates  are  quite  small 
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Fig.  1.  1-D  and  2-D  model  geometries  (not  to  scale). 


(25  ml  min-1)  and  allow  the  gases  to  come  up  to  cell  tem¬ 
perature  by  the  time  they  get  to  the  actual  cell.  The  above  was 
verified  by  modelling  only  the  fluid  flow  and  heat  transfer 
using  the  2-D  geometry. 

In  this  work,  the  one-dimensional  models  are  solved  using 
MATLAB  [20]  while  COMSOL  Multiphysics  [21]  is  used  to 
solve  the  two-dimensional  models.  MATLAB  is  a  general  pur¬ 
pose  computational  software  platform  and  language.  COMSOL 
Multiphysics  is  a  finite  element  method  based  modelling  pack¬ 
age  for  the  simulation  of  any  physical  process  that  can  be 
described  using  ordinary  or  partial  differential  equations.  It  has 
an  easy  to  use  graphic  user  interface  that  handles  all  modelling 
steps  from  geometry  generation  and  meshing,  to  defining  the 
differential  equations  and  boundary  conditions,  to  solving  and 
post-processing  the  results.  COMSOL  Multiphysics  has  a  MAT¬ 
LAB  interface  that  allows  easy  transfer  of  data  from  one  to 
the  other.  It  is  this  interface  that  allows  the  use  of  optimiza¬ 
tion  routines  available  in  MATLAB  for  the  parameter  estimation 
described  later. 


in  turn  requires  additional  kinetic  parameter  values  for  each 
reaction,  and  as  discussed  earlier,  these  are  unavailable  for  the 
anode  used  in  the  experiments.  This  is  also  the  reason  why 
reaction  (1)  is  the  only  electrochemical  reaction  modelled  in 
this  study. 

•  Nemst  equation  requires  the  compositions  of  the  products 
(H2O  and  S2  for  reaction  (1))  as  well  as  the  reactants.  These 
were  not  measured  during  the  experiments  that  are  being  mod¬ 
elled  in  this  study. 

2.2.  1-D  explicit  model 

In  developing  this  model  we  followed  a  similar  approach  to 
that  detailed  in  Kim  et  al.  [23]  and  Chan  et  al.  [24].  A  one¬ 
dimensional  analysis  was  done  along  the  thickness  of  the  cell 
(Fig.  1)  to  obtain  a  model  Eq.  (3)  that  relates  the  voltage  of  the 
working  cell,  Vcen  to  the  current  density,  i  drawn  from  the  cell: 

^cell  =  £oc  —  IRq  —  0?LCt,  anode  —  ^?act— cathode 

—  *7conc,  anode  —  *7  cone— cathode  (3) 


2.7.  Open  circuit  voltage 

In  most  fuel  cell  models,  the  open  circuit  voltage  or  the  equi¬ 
librium  voltage  across  a  fuel  cell  is  calculated  using  the  Nernst 
equation  [22].  In  all  the  models  discussed  in  this  work,  however, 
the  open  circuit  voltage,  Eoc  is  taken  from  experimental  data 
instead  of  being  calculated  using  the  Nernst  equation.  This  is 
done  for  the  following  reasons: 

•  The  open  circuit  voltage  in  a  H2S  SOFC  is  probably  a  mixed 
voltage  from  the  different  possible  reactions  at  the  anode 
[4,7,8].  Modelling  a  mixed  voltage  requires  that  the  kinetics 
of  each  possible  electrochemical  reaction  be  considered.  This 


where  Rq  is  the  total  ohmic  resistance  (cell  area  specific)  of  the 
cell,  rj act  the  activation  voltage  losses  at  the  anode  and  cathode, 
and  Oconc  are  the  concentration  (mass  transfer)  voltage  losses. 
The  activation  and  concentration  (mass  transfer)  voltage  losses 
in  our  model  are  given  in  Eqs.  (4)-(8)  below: 
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(7) 


in  the  anode,  cathode,  and  electrolyte.  The  model  domain  is  a 
line  that  extends  from  the  current  collector  on  the  anode  to  the 
current  collector  on  the  cathode  along  the  axis  of  the  cell  and 
includes  the  electrodes  and  the  electrolyte  (see  bottom  of  Fig.  1). 


RT  /  i  \ 

^conc— cathode  —  TT  In  1  1  —  -7  1 

4/7  V  lc,  l) 

(8) 

7*142  S,  in  4EDH2S,eff 

la’l~  (2+Pn2SM/P)  RTla 

(9) 

P02M  4FDo2.e„ 

(10) 

lc'1  -  (1  -  Po2,in/P )  RTlc 

where  l  is  the  length,  a  the  electrical  conductivity  of  the  elec¬ 
trode  or  electrolyte,  ne  the  number  of  electrons  transferred  in 
the  electrode  reaction,  and  P  is  the  total  or  absolute  pressure. 
For  both  the  1-D  models,  we  assume  that  P  is  constant  at 
1.013  x  105  Pa  in  both  electrodes.  The  ohmic  resistance,  given 
by  Eq.  (4),  includes  the  area  specific  resistance  of  the  zirconia 
electrolyte,  the  electrodes  and  the  contact  resistance.  Any  current 
lead  wire  or  inter-phase  resistance  is  bundled  into  /Contact-  The 
activation  losses  assume  a  Butler-Volmer  [22]  kinetic  expres¬ 
sion  where  the  charge  transfer  coefficients  are  0.5,  and  i®, 
are  the  exchange  current  densities.  /Th2o  in  Eq.  (7)  is  a  func¬ 
tion  of  initial  partial  pressures  of  H2O  and  H2S,  while  Ks2  is 
a  function  of  initial  partial  pressures  of  S2  and  F^S.  The  con¬ 
centration  voltage  loss  terms  have  the  limiting  current  densities 
ia,  1  and  ic,\  as  parameters,  which  in  turn  are  functions  of,  among 
other  things,  the  effective  diffusivities  of  the  reactants,  Z\e ff- 
Effective  diffusivity  is  a  function  of  the  bulk  diffusivity 
Knudsen  diffusivity  D^k,  and  the  ratio  of  porosity  to  tortuousity 
dr. 


Afc.eff 
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-1 
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The  solutions  to  the  mass  transfer  equations  for  the  reactants 
and  products  are  used  to  obtain  the  above  terms  for  the  con¬ 
centration  losses.  The  form  of  the  mass  transfer  equations,  and 
thus  the  solutions  are  identical  to  those  for  the  1-D  implicit 
model  discussed  in  the  next  section.  These  solutions  or  concen¬ 
tration  profiles  Eqs.  ((24)-(26)),  however,  need  to  be  linearized 
in  order  to  get  the  explicit  form  for  the  concentration  losses.  Of 
course,  the  1-D  explicit  model  does  not  consider  the  coupling 
between  the  activation  and  concentration  losses  as  the  form  of 
the  Butler-Volmer  equation  used  Eqs.  (5)  and  (6)  does  not  con¬ 
tain  any  concentration  terms,  while  the  1-D  implicit  model  does 
so  through  the  concentration  terms  in  the  electrochemical  rate 
equations  Eqs.  (20)  and  (21). 


2.3.  1-D  implicit  model 

The  governing  ordinary  differential  equations  (ODEs)  used 
in  our  1-D  implicit  model  are  presented  in  this  section.  A  1-D 
model  for  this  geometry  cannot  consider  flow  and  mass  transfer 
in  the  channels.  Thus,  the  model  equations  only  need  to  describe 
the  mass  transfer  in  the  electrodes  and  the  voltage  distribution 


2.3.1.  Governing  equations:  mass  transfer  and  voltage 

2. 3. 1.1.  Mass  transfer.  The  reactants,  oxygen  on  the  cathode 
side  and  hydrogen  sulfide  on  the  anode  side,  are  transported  from 
the  flow  channels  to  their  respective  electrode-electrolyte  inter¬ 
faces  where  the  reaction  occurs.  Similarly,  the  products  on  the 
anode  side  are  transported  from  the  anode-electrolyte  interface 
into  the  fuel  flow  channel.  This  mass  transfer  of  the  reactants  and 
products  is  modelled  using  the  convection-diffusion  equation: 

^  +  =0  (12) 


v  =  — 


& 

nk 


(13) 


where  Q  is  the  molar  concentration  of  the  kth  component,  Ct 
the  total  molar  concentration,  and  v  is  the  mole  average  velocity 
of  the  gas.  £={El2S,  H2O,  S2}  in  the  anode  and  {O2}  in  the 
cathode.  As  no  carrier  gas  is  used  on  the  anode  side  in  the  exper¬ 
iments,  we  can  model  transport  of  H2S  and  H2O  on  the  anode 
side  and  obtain  the  concentration  of  S2  at  any  point  by  using 
the  identity:  Cs2  =  Ct  —  Cr2s  —  Ch2o-  We  °btain  Eq.  (13)  for 
v  by  considering  the  overall  molar  flux  towards  or  away  from 
the  reaction  interface  and  assuming  that  the  total  pressure,  and 
thus  the  concentration  remains  constant  in  the  electrodes. 


2. 3. 1.2.  Potential  distribution.  The  fuel  reacts  on  the  anode- 
electrolyte  interface  and  releases  electrons.  These  electrons 
travel  through  the  anode  to  the  outer  circuit  and  come  around 
through  the  cathode  to  the  cathode-electrolyte  interface  where 
they  are  consumed  in  the  reaction  with  oxygen.  To  complete 
the  circuit,  the  oxide  ions  produced  at  the  cathode-electrolyte 
interface  travel  across  the  electrolyte  to  the  anode-electrolyte 
interface: 

H2S  +  O2-  &  H2O  +  ±S2  +  2e-  (14) 

02  +  4e“  <s>  202-  (15) 


The  voltage  and  current  distribution  in  the  electrodes  and  the 
electrolyte  due  to  this  electronic  and  ionic  transport  is  modelled 
using  Ohm’s  law: 


(16) 


where  oj  is  the  electrical  conductivity  and  <pj  is  the  electrical 
potential  of  the /h  sub-domain  (anode,  electrolyte,  or  cathode). 


2.3.2.  Boundary  conditions 

To  get  the  concentration  and  voltage  profiles,  the  following 
boundary  conditions  need  to  be  solved: 
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•  Concentrations  of  O2  in  air,  and  H2S,  H2O  in  fuel  given  at 
the  inlets: 

CH2sla^fUei,iniet  =  ^H2S,in 

^H20  1 9 £2 fuel, inlet  =  CH20,in  (17) 

Q)2laftair,inlet  =  ^02,in 

•  Molar  flux,  of  the  reactants  into,  and  the  products  out  of, 
the  electrolyte  governed  by  the  local  current  density  at  the 
electrode-electrolyte  interface: 


—  £>H2S,eff 

—  £>H20,eff 
“A)2,eff 
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d  z 
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d  z 
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9  £2  anode ,  electrolyte 


9  £2  anode ,  electrolyte 


3^c  athode ,  elec  troly  te 


1  anode 

2  F 

1  anode 

'  2F~ 

1  cathode 

4F 

(18) 


•  The  voltage  at  the  anode  and  cathode  current  collectors  is 
specified: 


d  (pa 

1  1 

dz 

d z 


=  (<Pa  -  0) 


9  ^  anode ,  collector 


^  £2  cathode ,  collector 
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111 


(19) 


where  [all]  contact  is  the  reciprocal  of  the  area  specific  contact 
and  lead  resistance  for  each  electrode. 

•  Current  density  at  the  electrode-electrolyte  interfaces  is  given 
by  the  electrochemical  reaction  rate,  which  is  dependent  on 
the  local  concentrations  of  the  active  species  and  the  local 
voltage  difference  between  the  electrode  and  the  electrolyte 
[22]: 


electrochemical  rate  equation  depends  on  local  species  concen¬ 
trations  as  well  as  the  local  potential  difference  between  the 
electrode  and  electrolyte  phases. 

Solution  of  the  mass  transfer  ODEs  and  boundary  conditions 
gives  the  component  partial  densities  and  electric  potentials  as 
functions  of  current  density  and  distance  from  the  anode  current 
collector  along  the  cell  axis: 


Ch2S  =  “2 Cf  +  (2 Cf  +  CH2S,in)  exp 


iz 


4CVF£>H2S,eff 


(24) 


Cr2o  =  2 C*  +  (-2 C*  +  CH2o,in)  exp 


iz 


4QF£>H2o,eff 


(25) 


Co2  =Ct  +  (-Cf  +  Co2,in)  exp 


i\z  —  (l a  +  Ijn  +  4)]  \ 

4QFZ)02,eff  / 

(26) 


4>j  =  nijZ  +  Cj  (21) 

where  mj  and  cj  are  integration  constants  for  the  voltage  profile 
equation  for  the  jth  sub-domain  (anode,  electrolyte,  or  cathode). 

Although  we  obtain  the  concentration  profiles  above,  they 
are  functions  of  the  current  density  which  in  turn  is  a  function 
of  the  concentrations  at  the  inter  phase  boundaries.  The  elec¬ 
tric  potential  boundary  conditions  Eqs.  ((19)-(23))  provide  a 
set  of  coupled  non-linear  algebraic  equations  that  cannot  be 
solved  analytically.  In  this  work,  these  six  equations  in  six 
unknowns  (my  and  cy),  are  solved  numerically  using  MATLAB 
[20]. 


-<* a 
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0  c  —  0cathode  —  0electrolyte  —  ^0c  (22) 

where  and  are  the  exchange  current  densities  of  the  anode 
and  cathode,  /3  the  charge  transfer  coefficient,  Q,ref  the  reference 
(in  this  case,  inlet)  component  densities  and  A0^  and  A0*?  are 
the  anode  and  cathode  equilibrium  potentials. 

The  boundary  conditions  given  by  Eqs.  (20)  and  (21)  are 
called  electrochemical  rate  equations  because  while  a  chemi¬ 
cal  rate  equation  depends  on  local  species  concentrations,  an 


2.4.  2-D  models 

The  first  2-D  model,  built  using  COMSOL  Multiphysics 
[21],  solves  the  (i)  non-isothermal  flow  equations  and  Darcy’s 
law  for  velocity  and  pressure  in  the  flow  channels  and  elec¬ 
trodes,  respectively,  (ii)  convection-diffusion  mass  transfer 
equations  for  partial  densities  of  H2S,  H2O,  S2  on  the  fuel 
side  and  O2  on  the  air  side  and,  (iii)  Laplace  equation  (Ohm’s 
law)  for  the  voltage/current  distribution  in  the  electrodes  and 
electrolyte. 

The  second  2-D  model  uses  the  Maxwell-Stefan  equations 
(instead  of  Fick’s  law)  to  model  the  multicomponent  mass  trans¬ 
fer  in  the  fuel  cell. 
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2.4.1.  Flow  in  gas  channels 

The  fluid  flow  in  the  fuel  and  air  channels  is  important  as 
it  describes  the  transport  of  reactants  to  and  the  products  from 
the  cell.  The  non-isothermal  flow  mode  is  used  here  because 
the  composition  changes  in  the  fuel  and  air  channels  lead  to 
small  density  variations.  These  density  variations  are  taken  into 
account  by  adding  a  term  to  the  momentum  equations  and  mod¬ 
ifying  the  continuity  equation  as  shown  in  the  governing  PDEs 
for  the  flow  fields: 


-V  • 


MVv  +  (Vv)T) 


|/*(V  •  v)I 


+  p(v  •  Vv)  +  Vp  =  0, 


V  •  (pv)  =  0 


(28) 


where  /x  is  the  viscosity  of  the  fluid,  v  the  velocity  vector,  p  the 
density,  p  the  pressure,  and  I  is  the  identity  matrix. 


2.4. 1.1.  Boundary  conditions  (channel  flow). 


•  Zero  velocity  (no-slip)  at  the  walls: 


2.4.2. 1.  Boundary  conditions  ( electrode  flow ). 


•  Pressure  at  each  electrode-flow  channel  interface  is  equal  to 
the  pressure  in  the  flow  channel  at  that  interface: 


P anode  I  3  £2 anode, fuel  channel  Pfuel  channel  I  3^ anode, fuel  channel 

Pcathode  I  3£2cathode,  air  channel  Pair  channel  I  3^cathode,  air  channel 


(35) 

(36) 


•  Flow  into  the  electrodes  at  the  electrode-electrolyte  inter¬ 
face  is  the  net  difference  in  the  amount  of  the  products  being 
formed  and  the  reactants  being  consumed  at  that  interface: 


-n.v|an, 


'electrode,  electrolyte 


1  /  /electrode 

p  ^  F 


%k  M  K  .  \ 

?7 

(37) 


where  n  is  the  normal  vector,  /electrode  the  current  density,  F 
the  Faraday’s  constant,  the  stoichiometric  coefficient,  M k 
the  molar  weight,  %  the  number  of  electrons  transferred  per 
molecule,  and  is  the  diffusive  flux  of  the  kth  component. 

•  Radial  symmetry  along  the  axis  of  the  cell: 


^  t)  1 3£2wans 


(29) 


•  Velocity  continuous  across  the  channel-electrode  interfaces: 


vflow  channels 


^electrodes 


3 £2 flow  channels— electrodes 


(30) 


•  Fully  developed  laminar  flow  at  the  fuel  and  air  inlets.  Given 
the  volumetric  flow-rate,  the  axial  velocity  profile  is  calcu¬ 
lated  as: 


u  —  0, 


(31) 


where  u  and  v  are  the  r  and  z  components  of  the  velocity 
vector,  V  the  temperature  corrected  volumetric  flow-rate,  and 
ri  is  the  inner  radius  for  the  inlet  tube. 

•  Pressure  specified  and  flow  normal  to  the  boundary  at  the 
outlets  of  the  flow  channels: 


p  =  0,  t  •  v  =  0  (32) 

where  t  is  the  tangential  vector. 

•  Radial  symmetry  along  the  axis  of  the  flow  channels: 


dv 

dr 


=  0, 

r=0 


U  \r=0  —  0 


(33) 


2.4.2.  Flow  in  electrodes 

The  flow  in  the  porous  electrodes  is  modelled  using  Darcy’s 
law: 

v  =  —  —  Vp,  V  •  (pv)  =  0  (34) 

li 


dp 

dr 


=  0 


r= 0 


(38) 


2.4.3.  Mass  transfer 

Two  mass  transfer  models  are  used  for  the  2-D  fuel  cell  mod¬ 
els  in  this  work.  The  convection-diffusion  Eq.  (39)  is  used  in 
both  models  to  descibe  the  mass  transfer  of  the  reactants  and 
products  in  the  flow  channels  and  the  electrodes: 

V  •  jf  +  pv  •  Vtu;  =  0  (39) 

The  diffusive  flux  is  given  by  Fick’s  law  Eq.  (40)  in  the  first 
2-D  model,  and  by  the  Maxwell-Stefan  diffusive  flux  Eq.  (41) 
in  the  second  [25]: 

j.  =  -pDiV(wi)  (40) 

where  j/  is  the  diffusive  flux  of  /,  p  the  density  of  the  gas  mix¬ 
ture  and  Wi  is  the  mass  fraction  of  the  /th  component.  D;  is  the 
diffusivity  of  component  /:  bulk  diffusivity  in  the  flow  channels, 
effective  diffusivity  in  the  electrodes.  Mass  fractions  are  used 
instead  of  mole  fractions  or  molar  concentrations  in  the  Fick 
diffusive  flux  term  because  the  mass  transfer  is  coupled  to  the 
fluid  flow  equations  where  the  fluid  velocity  v  is  a  mass  average 
quantity: 


In  Eq.  (41),  n  is  the  total  number  of  species  in  the  mixture,  Xj 
is  the  mole  fraction  of  species  j,  and  Dp  is  the  multicomponent 
composition  dependent  diffusivity  of  species  i  in  j,  which  is 
given  by: 

Xi%k  ^i) jk 

-  =  —'WiWk„  ~ - , 

D,k  k  V ( adj  Bi  )jk 

(Bi)kj  =  Dkj  —  Di j ,  i  j 


where  k  is  the  permeability  of  the  porous  electrode. 


(42) 
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where  Du  are  multicomponent  Maxwell-Stefan  diffusivities 
that  for  low  density  gas  mixtures  can  be  approximated  by  com¬ 
position  independent  binary  diffusivities,  and  can  be  estimated 
using  the  Fuller,  Schettler  and  Giddings  equation.  The  values 
used  in  this  work  are  given  in  Section  2.5. 

Mass  fractions  and  mole  fractions  are  related  through: 


XjMj 


(43) 


For  Fick’s  law,  effective  diffusivity  in  the  electrodes  is  defined 
in  Eq.  (11),  but  for  Maxwell-Stefan  mass  transfer,  Eq.  (44)  is 
used: 


2.4.4.  Voltage  and  current  distribution 

The  voltage  and  current  distribution  in  the  electrodes  and  the 
electrolyte  due  to  this  electronic  and  ionic  transport  is  modelled 
using  the  vector  form  of  Ohm’s  law: 


^  ’  (  °rm^70m)  —  0 

v  .  (-crflV0fl)  =  0  (49) 

V  •  (— crcV0c)  =  0 

where  crm,  a a,  oc  are  electrical  conductivities  of  the  electrolyte, 
anode  and  cathode,  and  0m,  <j>a,  <pc  are  electrical  potentials  of 
the  electrolyte,  anode,  and  cathode. 


Dik,e&=-Dik  (44) 

r 

As  discussed  earlier,  no  carrier  gas  is  used  on  the  anode  side  in 
the  experiments.  Thus,  we  can  model  transport  of  F^S  and  F^O 
on  the  anode  side  and  obtain  the  mass  fraction  of  S 2  at  any  point 
by  using  the  identity:  wS2  =  1  —  wu2s  —  wh2o 


2.4.3. 1.  Boundary  conditions  (convection-diffusion  equation). 

•  Mass  fraction  of  O2  in  air,  and  H2S,  S2,  H2O  in  fuel  given  at 
the  flow  channel  inlets: 

WH2sl3S2fueUnie,  =  WH2S,in 

wH20  1 3£2fuei, inlet  =  wH20,in  (45) 

W02l3Qair,iniet  =  Wo2,in 

•  Zero  flux  at  the  inner  and  outer  tube  walls: 


2. 4.4.1.  Boundary  conditions  ( Ohm ’s  law ). 


•  The  voltage  at  the  anode  and  cathode  current  collectors  is 
specified: 


n  •  ianode  1 3^ anode  collector  (0a  -  0) 

^  ^cathode  1 3£2cathode,  collector  (*Pc 
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contact 
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(50) 


where  i  is  the  current  density  vector  [cr//]contact  is  the  area 
specific  contact  and  lead  conductance. 

•  Electrical  insulation  at  the  outer  tube  walls: 


n  '  ianode  I  Stalls  =  ^ 

n  ’  ^cathode  1 3£2Walls  =  0  (51) 

n  '  ^electrolyte  I  Stalls  —  0 


n  •  N^|a^walls  =  0 
=  j*  +  pwk\ 


(46) 


where  N&  is  the  total  mass  flux  of  the  kth  component. 


•  Current  density  at  the  electrode-electrolyte  interfaces  is  given 
by  the  electrochemical  reaction  rate,  which  is  dependent  on 
the  local  concentrations  of  the  active  species  and  the  local 
voltage  difference  between  the  electrode  and  the  electrolyte: 
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(53) 


•  Mass  flux,  of  the  reactants  out  of,  and  the  products  into, 
the  electrodes  governed  by  the  local  current  density  at  the 
electrode-electrolyte  interface: 


®  NH2  S  1 3£2  ano(ie  5  electrolyte 
n  •  NHzOls^anode  ,  electrolyte 
^  NO,  1 3^cathode  ,  electrolyte 


^H2S*  anode 

IF 
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2  F 
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4 F 


(47) 


da  —  0anode  ^electrolyte 

0c  —  0cathode  —  0electrolyte  —  A0C  (54) 

where  and  are  the  exchange  current  densities  of  the  anode 
and  cathode,  f  the  charge  transfer  coefficient,  xk  the  mole 
fractions  at  the  boundary,  x^ef  the  reference  (in  this  case, 
inlet)  mole  fractions  and  A 0^  and  A  0^  are  the  anode  and 
cathode  equilibrium  potentials. 

•  Radial  symmetry  along  the  axis  of  the  electrodes  and  the  elec¬ 
trolyte: 


•  Radial  symmetry  along  the  axis  of  the  flow  channels  and  the 
electrodes: 


dwk 

dr 


=  0 

r=0 


(48) 


90/ 

dr 


=  0 

r=0 


(55) 


The  mesh  generated  for  the  2-D  models  above  had  3853  trian¬ 
gular  elements,  which  translates  to  30,229  degrees  of  freedom. 
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We  used  the  non-linear  parametric  solver  (with  the  UMFPACK 
direct  solver  for  the  linear  subsystem)  in  COMSOL  Multiphysics 
to  solve  the  model  at  equally  spaced  voltage  intervals  from  near 
OCV  to  0  V.  Grid  independence  of  the  results  was  insured  by 
comparing  the  model  output  for  the  current  density  at  0  V  and 
850  °C,  against  a  mesh  4  times  denser.  The  relative  error  in  the 
calculated  current  denstiy  was  of  the  order  of  10-6  between  the 
mesh  used  in  this  work  (with  3853  elements)  and  the  test  mesh 
with  15,412  elements. 

2.5.  Parameters  used  in  the  models 

There  are  a  number  of  parameters  in  the  models  presented 
above  and  these  parameters  need  to  be  estimated  before  the  mod¬ 
els  can  be  solved.  Physical  constants,  such  as  the  gas  constant 
R ,  and  the  values  assigned  to  them  in  the  models  are  summa¬ 
rized  in  Nomenclature.  Transport  parameters  needed  include 
diffusivities,  viscosities,  and  electrical  conductivities,  while  the 
electrochemical  parameters  needed  are  equilibrium  electrode 
potential,  exchange  current  density  and  the  charge  transfer  coef¬ 
ficient  for  each  electrode. 

Binary  diffusion  coefficients  Dy  are  estimated  using  the 
Fuller,  Schettler  and  Giddings  relation  given  in  Reid  et  al.  [26]. 
Knudsen  diffusion  coefficients  Dk,Jc  are  calculated  assuming 
smooth  round  pores.  Viscosity  of  the  gases  on  both  sides  of 
the  fuel  cell  /xa ir,/Xfuei  is  also  estimated  using  Brokaw’s  method 
outlined  in  Ref.  [26] .  The  densities  of  the  gases  are  calculated 
using  the  ideal  gas  law  and  the  mass  fractions  of  the  compo¬ 
nents.  Experimental  values  for  electrode  permeability  k  were 
not  available  and  order  of  magnitude  estimates  are  used  based 
on  permeability  values  calculated  for  close  packings  of  spheres 
[27]. 

Electrical  conductivity  of  the  anode  materials  oa  used  was 
measured  by  Liu  [28]  and  those  values  are  used  in  this  work. 
The  cathode  used  in  the  experiments  was  Pt,  and  the  value  for 
Pt  conductivity  at  the  operating  temperatures  is  corrected  for 
porosity  using  relations  given  in  Ref.  [29].  Values  for  conductiv¬ 
ity  of  the  YSZ  electrolyte  am  are  obtained  using  the  temperature 
dependent  form  given  in  Ref.  [30]. 

The  equilibrium  potential  of  the  cathode  A  0^  is  set  to  the 
value  of  the  open  circuit  potential  at  each  operating  tempera¬ 
ture  while  the  equilibrium  potential  of  the  anode  is  set  to 
zero.  The  open  circuit  potential  at  each  operating  temperature  is 
taken  from  the  experimental  i-V  data  sets.  The  charge  transfer 
symmetry  coefficients  ft  in  the  electrochemical  rate  equations 
Eqs.  (20),  (21),  (52)  and  (53)  are  assigned  the  default  value  of 
0.5  [22].  As  experimental  values  for  exchange  current  densities 
are  not  available,  we  estimate  them  using  our  models  and  exper¬ 
imental  i-V  data.  The  estimation  procedure  is  explained  in  the 
next  section. 

The  overall  wire  lead  and  contact  resistance  is  calculated 
by  comparing  the  cell’s  i-V  curve  to  the  IR  compensated  i-V 
curve  at  the  different  operating  temperatures.  The  overall  cell 
area  specific  resistance  is  given  by  Eq.  (56),  where  Vir  is  the  IR 
compensated  cell  voltage  and  V  is  the  actual  cell  voltage  when 
a  current  density  of  i  is  drawn  from  the  cell.  Four  such  values 
taken  at  different  points  of  the  i-V  curve  are  averaged  and  the 


Table  1 


Temperature  independent  parameters  in  the  models 


Parameter 

Value 

P  (Pa) 

1.013  x  105 

Vfuel  (ml  min-1) 

25 

Vair  (ml  min-1) 

25 

*H2S,in 

0.985 

*H20,in 

0.01 

■^02,in 

0.21 

e 

0.4 

X 

4 

r pore,av  (m) 

2  x  10-6 

la  (m) 

10-4 

lc  (m) 

10-4 

lm  (m) 

2  x  10-4 

Ka  (m2) 

10-14 

Kc  (m2) 

10-14 

known  resistances  of  the  different  cell  layers  are  subtracted  from 
the  overall  resistance  to  obtain  the  overall  wire  lead  and  contact 
resistance  for  the  cell  Eq.  (57).  This  overall  resistance  is  then 
divided  into  equal  halves  in  our  models  and  placed  at  the  current 
collectors  for  each  electrode: 


Vir  —  V 

Rq  =  —. - 

i 


(56) 


l 


o 


contact,  t 


=  Rn- 


~r 

~r 

cr 

YSZ 

< 7 

anode 

a 

cathode 

(57) 


/ 

.  **  J  contact 


/ 

J  contact,? 


(58) 


The  numerical  values  of  all  the  parameters  used  are  given  in 
Tables  1  and  2. 


3.  Parameter  estimation 

As  discussed  earlier,  the  electrochemical  parameters  on  the 
anode  side  are  completely  unknown  whereas  those  on  the  cath¬ 
ode  side  are  uncertain.  In  this  work  we  estimate  exchange  current 
densities  and  using  non-linear  least  squares  [31].  This 


Table  2 

Temperature  dependent  parameter  values  used 


Parameter 

750  °C 

800  °C 

850  °C 

£>h2s,h2o  (m2) 

1.89  x  10“4 

2.05  x  10“4 

2.23  x  10“4 

£>h2s,s2  (m2) 

9.12  x  10-5 

9.91  x  10“5 

1.07  x  10“4 

£>h2o,s2  (m2) 

1.42  x  10“4 

1.54  x  10“4 

1.67  x  10“4 

£>o2,n2  (m2) 

1.8  x  nr4 

1.96  x  10“4 

2.12  x  10-4 

Dk,  h2s  (m2) 

1.06  x  10“3 

1.09  x  10“3 

1.12  x  10“3 

Dk,  h2o  (m2) 

1.46  x  10“3 

1.50  x  10“3 

1.53  x  10“3 

Dk, s2  (m2) 

7.76  x  10“4 

7.94  x  10“4 

8.13  x  10“4 

Dk,  o2  (m2) 

1.1  x  10“3 

1.12  x  10“3 

1.15  x  10“3 

/Zair  (kg  (ms)-1) 

3.71  x  10“5 

3.82  x  10“5 

3.93  x  10“5 

Mfuel  (kg  (rns)-1) 

3.7  x  10“5 

3.85  x  10“5 

4.01  x  10“5 

A  4>°c  (V) 

0.886 

0.880 

0.872 

cra  ((£2  m)-1) 

13.1 

9.88 

5.26 

crc  ((£2  m)-1) 

1.11  x  106 

1.06  x  106 

1.02  x  106 

crm  ((£2m)-1) 

1.49 

2.39 

3.66 

[cr//] contact  (£2  1  m  2) 

7670 

7000 

7730 
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method  minimizes  the  sum  of  the  squares  of  the  difference 
between  experimental  data  points  and  the  model  output  by  vary¬ 
ing  the  parameters  to  be  estimated. 

The  data  are  filtered  before  being  sent  to  the  estimation  algo¬ 
rithm.  This  filtering  includes: 

(1)  The  i-V  data  are  smoothed  to  remove  noise.  The  Robust 
Lowess  smoothing  method  provided  in  MATLAB’s  Curve 
Fitting  Toolbox  was  used  [32]. 

(2)  Repeating  data  points  and  data  values  against  the  trend  are 
removed. 

(3)  Current  density  values  are  interpolated  for  an  equally  spaced 
vector  of  voltage  values  to  get  a  i-V  matrix  of  a  manageable 
size.  For  all  four  models  used  in  this  study,  the  V  vector  was 
defined  by  16  equally  spaced  points  that  went  from  99%  of 
the  open  circuit  potential  of  the  cell  to  0  V. 


3.3.  Computational  resources  used 

All  parameter  estimation  routines  were  run  using  the  function 
Isqnonlin  in  the  MATLAB  Optimization  Toolbox  [31].  The  1- 
D  explicit  model’s  parameter  estimation  took  ~1  s  on  a  3  GHz 
Intel  Pentium  computer.  The  1-D  implicit  model’s  parameter 
estimation  took  ~5  min  while  parameter  estimation  for  the  2-D 
models  took  from  5  to  14  h  per  data  set.  As  seen  in  the  section 
on  modelling,  the  1-D  explicit  model  calculates  the  output  Vm(i) 
as  an  explicit  function  of  i,  while  the  1-D  implicit  model  needs 
to  solve  a  system  of  6  coupled  non-linear  algebraic  equations  to 
calculate  im(V).  For  the  2-D  models  presented  here,  COMSOL 
Multiphysics  solves  for  ~30,000  degrees  of  freedom  to  calculate 
the  current  drawn  by  the  cell  at  a  given  cell  voltage. 

4.  Results  and  discussion 


The  parameter  estimation  routine  takes  this  matrix  of  i-V  data 
at  each  temperature  and  uses  the  models  to  iterate  to  z'o  values 
that  fit  the  given  data  best.  The  1-D  explicit  model  calculates  the 
cell  voltage  as  a  function  of  specified  average  current  density 
while  the  1-D  implicit  and  the  2-D  models  calculate  current 
density  for  a  specified  cell  operating  voltage.  Thus,  the  parameter 
estimation  needs  to  be  set  up  differently  depending  on  the  model 
being  fitted. 


3. 1 .  Mathematical  formulation  for  1  -D  explicit  model 


n 

min yy  Vm(/*,  T)  -  14. r I2 

iafc  fc=l 

subject  to  : 

Vm  =  f(i,  T,  Z°,  fic,  P) 


(59) 


where  Vm( i^T)  is  the  output  cell  voltage  given  by  the  model  for 
the  current  density  4  while  Vkj  is  the  experimentally  observed 
cell  voltage  at  current  density  4-  The  vector  p  represents  the 
known  parameters  (e.g.,  oj,  Dj  e ff)  needed  to  solve  the  model  for 
Vm,  and  n  is  the  number  of  data  points. 


3.2.  Mathematical  formulation  for  1-D  implicit  and  2-D 
models 


n 

min ^2{im(Vk,  T)  -  ikj}2 

iaJck=l 

subject  to  : 

im  =  g(V,  T,  i%  i°c,  p*) 


(60) 


where  im(Vk,T)  is  the  output  cell  current  density  given  by  the 
model  for  a  cell  voltage  of  14  and  4,r  is  the  experimentally 
observed  current  density  for  the  same  cell  voltage  of  14 .  The 
vector  p*,  again  represents  the  known  parameters  needed  to 
solve  the  model  for  im. 


The  results  of  the  parameter  estimation  for  all  four  mod¬ 
els  are  summarized  in  Table  3.  The  last  column  on  the  right 
hand  side  gives  the  scaled  fit  for  all  the  models  at  the  three 
temperatures.  These  numbers  are  obtained  by  scaling  the  nor¬ 
malized  fit  with  the  maximum  value  in  the  data  being  fitted: 

—  14)2/n)/(max(14))  for  the  1-D  explicit  model  and 

\J (52  dm  —  4)2/w)/(max(4))  f°r  the  other  two  models. 

There  are  several  trends  that  can  be  discerned  from  this  data. 
The  first  one  is  that  i®  and  i®  increase  with  increasing  temper¬ 
ature  for  all  the  models.  This  trend  is  expected  for  catalyzed 
electrochemical  reactions. 

Model  fit  to  data  improves  with  increasing  temperature  for 
all  the  models.  The  models  fit  the  data  at  850  °C  better  than  at 
800  and  750  °C  (see  Figs.  2-4).  A  simplistic  explanation  for 
this  trend  is  that  the  performance  curves  become  less  non-linear 
and  thus  easier  to  fit,  with  increasing  temperature.  However, 
the  models  do  allow  a  non-linear  i-V  response  through  the 
electrochemical  rate  equations  Eqs.  (52)  and  (53).  We  believe 
this  mismatch,  more  severe  at  lower  temperatures  due  to  the 
lower  i°j,  is  related  to  the  symmetric  charge  transfer  coefficients 
(2(1  —  f)  and  2/3  with  f>  =  0.5)  used  in  the  electrochemical  rate 
equations.  These  values  are  obtained  by  assuming  a  reaction 


Table  3 

Parameter  estimation  results 


Model  used 

r(°c) 

*2  (Am”2) 

(Am“2) 

Scaled  fit 

1-D  explicit 

750 

70.0 

70.0 

2.85  x  10"2 

800 

180.7 

180.6 

1.50  x  10"2 

850 

344.9 

344.9 

3.98  x  10“3 

1-D  implicit 

750 

60.1 

62.3 

4.14  x  10“2 

800 

157.5 

163.3 

2.15  x  10"2 

850 

300.0 

309.1 

9.32  x  10“3 

2-D  Fickian  mass 

750 

55.0 

73.3 

4.14  x  10"2 

transfer 

800 

181.7 

187.2 

2.51  x  10“2 

850 

571.1 

356.3 

1.59  x  10"2 

2-D  Maxwell-Stefan 

750 

54.1 

74.5 

4.13  x  10"2 

mass  transfer 

800 

163.9 

199.3 

2.50  x  10“2 

850 

565.8 

362.9 

1.60  x  10"2 
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Model  output  using  estimated  i°  at  1023  K 


Fig.  2.  Model  output  using  optimized  z°  at  750  °C. 


mechanism  with  a  rate  controlling  two  electron  transfer  step. 
Ongoing  research  in  our  group  is  focussing  on  deriving  better 
electrochemical  rate  equations. 

According  to  the  scaled  fit  values  in  Table  3  the  1-D  explicit 
model  seems  to  give  a  better  fit  to  the  data  at  all  three  tempera¬ 
tures.  It  does  seem  counterintuitive  that  a  model  that  we  know  to 
be  too  simple  to  explain  all  that  is  happening  in  the  fuel  cell  gives 
a  better  fit  to  the  experimental  data.  However,  visual  examina¬ 
tion  of  the  model  fit  for  the  different  models  to  the  experimental 
i-V  data  does  not  show  any  significant  difference  between  the 
fit  for  the  1-D  versus  2-D  models,  and  as  we  note  later,  the  1- 
D  models  cannot  correctly  calculate  the  coupled  activation  and 
concentration  losses,  especially  at  the  low  flow-rates  used  here. 

All  models  give  values  for  i®  and  fic  that  give  a  fair  Arrhenius 
temperature  dependence.  The  Arrhenius  plot  for  the  and  fic 
given  by  the  2-D  model  using  Maxwell-Stefan  mass  transfer  is 


Model  output  using  estimated  i°  at  1073  K 


Model  output  using  estimated  i°  at  1 1 23  K 


Fig.  4.  Model  output  using  i°  at  850  °C. 


2-D  model  with  Maxwell-Stefan  mass  transfer 


Fig.  5.  Arrhenius  plot  for  i°  from  2-D  M-S  model. 


shown  in  Fig.  5.  This  plot  gives  the  fi  versus  1  IT  curves  obtained 
with  pre-exponentials  and  activation  energies  calculated  using 
the  optimal  and  fic.  The  plot  also  shows  the  actual  optimal 
values  for  the  exchange  current  densities  used  to  calculate  the 
least  squares  fitted  values  for  Aj  and  Ej  in  Eq.  (61).  These  pre¬ 
exponential  factors  and  activation  energies  are  given  in  Table  4: 

=  Aae~Ea/RT 


Table  4 

Arrhenius  parameters  estimated  for  i°  in  2-D  Maxwell-Stefan  model 


i° 

1° 

la 

lC 

A  (Am-2) 

1.41  x  1013 

4.35  x  109 

E  (kj  mole-1) 

224 

152 
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4.1.  The  1-D  models 

The  main  difference  between  the  1-D  explicit  and  the  1- 
D  implicit  models  is  that  in  the  explicit  model  the  reaction 
rate  is  independent  of  the  local  species  concentrations  (Eqs. 
(5)  and  (6))  whereas  the  reaction  rate  in  the  implicit  model 
depends  on  the  local  concentrations  at  the  reaction  surface  (Eqs. 
(20)  and  (21)).  In  the  implicit  model  the  reaction  rate  is  thus 
coupled  to  the  mass  transfer  whereas  the  activation  and  con¬ 
centration  losses  are  decoupled  and  separable  in  the  explicit 
model.  The  above  difference  between  the  two  models,  how¬ 
ever,  gives  no  clear  indication  as  to  how  it  would  affect  P 
estimates. 

4.2.  1-D  implicit  model  versus  2-D  Fick  model 

For  all  the  three  data  sets  or  temperatures,  the  P  values  in 
the  1-D  implicit  model  are  lower  than  those  in  the  2-D  model 
using  Fick  mass  transfer.  This  difference  can  be  explained  by 
the  difference  in  species  concentrations  at  the  interface  between 
the  electrodes  and  the  electrolyte  in  the  two  models.  The  2-D 
model  accounts  for  flow  and  mass  transfer  in  the  gas  channels 
and  is  thus  able  to  include  the  mass  transfer  resistance  in  the  flow 
domain.  The  1-D  implicit  model  cannot  account  for  the  inher¬ 
ently  two-dimensional  flow  field  in  a  button  cell  and  cannot 
account  for  mass  transfer  in  the  channels  correctly.  Therefore, 
the  parameter  estimation  for  the  1-D  implicit  model  ends  up 
assigning  a  higher  activation  resistance  (lower  P).  A  plot  of  the 
mass  fraction  profile  of  O2  along  the  axis  of  the  cell  is  given 
in  Fig.  6.  This  supports  the  above  discussion  by  illustrating  the 
difference  between  the  1-D  and  2-D  models  in  predicting  con¬ 
centrations  at  the  cathode-electrolyte  boundary  (surface  where 
the  cathode  reaction  occurs). 

The  above  discussion  also  implies  that  if  sufficiently  high 
flow-rates  are  used,  the  1-D  implicit  model  will  approximate  the 
2-D  model  better.  This  is  verified  in  Fig.  7  where  two  i-V curves 


02  mass  fraction  along  cell  axis 


Fig.  6.  O2  mass  fraction  profiles  along  axis  for  the  1-D  implicit  and  2-D  models 
at  850  °C. 


2-D  model  output  using  i°  from  fitted  1-D  model  at  850°C,  effects  of  varying  flowrate 


Fig.  7.  2-D  model  output  using  i°  from  1-D  model  at  850  °C  at  different  flow- 
rates. 


generated  using  the  2-D  Maxwell-Stefan  model  are  compared 
to  the  experimental  data  at  850° C.  The  two  curves  given  have 
the  same  values  for  Pa  and  Pc :  the  parameter  estimates  for  the  1-D 
implicit  model.  The  flow-rates  of  the  gas  streams,  however,  are 
25  ml  min-1  for  one  curve  and  250  ml  min-1  for  the  other.  As 
clearly  seen  in  Fig.  7,  the  model  output  for  the  higher  flow-rate  is 
able  to  approximate  the  experimental  data  quite  well  even  though 
the  kinetic  parameters  used  are  obtained  from  the  1-D  model. 
This  is  because  the  higher  flow-rate  dramatically  decreases  the 
mass  transfer  resistance  in  the  gas  channels. 

The  radial  variation  of  species  concentrations  is  insignificant 
as  seen  in  Fig.  8.  The  rise  in  wo2  from  the  left  to  the  right  of 
Fig.  8  is  due  to  diffusion  into  the  electrode  from  the  electrode 
edge  at  the  right. 


Fig.  8.  O2  mass  fraction  radial  profile  in  the  2-D  Maxwell-Stefan  model  at 
850  °C. 
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2-D  model  performance  curve  and  voltage  losses 


mass 

0.036 

0.034 

0.032 

0.03 

0.028 

0.026 

0.024 


fractions  of  H_2S  (top  half)  and  0_2  (bottom  half)  near  the  cell  assembly 


x  lO*^ 


Fig.  9.  Voltage  loss  profiles  in  the  2-D  Maxwell-Stefan  model  at  850  °C. 

4.3.  The  2-D  models 

The  2-D  models  used  differ  only  in  how  they  compute  the 
mass  transfer  of  the  different  species.  The  Maxwell-Stefan  mass 
transfer  formulation  is  superior  to  the  Fick  mass  transfer  form 
because  it  correctly  accounts  for  the  variation  in  the  diffusivites 
of  the  different  components  with  composition  [33]. 

Taking  the  2-D  model  with  Maxwell-Stefan  mass  transfer  to 
be  the  most  faithful  model,  we  now  discuss  the  profiles  of  the 
different  voltage  losses  in  the  SOFC.  Fig.  9  gives  the  cell  voltage 
as  well  as  various  voltage  losses  (commonly  called  overpoten¬ 
tials)  as  a  function  of  cell  current  density  according  to  the  second 
2-D  model  at  850  °C. 

The  plot  clearly  shows  that  the  highest  voltage  loss  is  due  to 
the  contact/current  lead  resistances  which  account  for  2.68  Q  or 
80%  of  the  total  ohmic  resistance  of  3 .33  Q  at  850  °C.  The  anodic 
and  cathodic  activation  losses,  respectively,  are  next  highest, 
followed  by  the  ohmic  loss  in  the  electrolyte.  The  ohmic  loss 
in  the  anode  is  roughly  14  times  smaller  than  the  ohmic  loss  in 
the  electrolyte,  while  that  in  the  cathode  is  another  4  orders  of 
magnitude  smaller  than  in  the  anode.  The  relative  magnitudes 
of  the  different  ohmic  losses  can  also  be  readily  calculated  from 
the  resistivity  values  and  the  breadths  of  the  different  phases 
(Tables  1  and  2). 

Fig.  10  gives  the  mass  fractions  of  the  reactants  near  the 
electrode-electrolyte  assembly  in  the  modelled  SOFC  at  an 
operating  voltage  of  0  V  (short-circuited  cell  drawing  maximum 
current).  The  vertical  boundary  on  the  left  is  the  axis  of  symme¬ 
try  of  the  cell,  the  gap  in  the  middle  is  the  electrolyte,  and  the 
top  half  is  the  fuel  side  while  the  bottom  half  the  air  side.  The 
mass  fraction  of  F^S  on  the  fuel  side  and  the  mass  fraction  of  O2 
on  the  air  side  is  shown  using  a  colour  scale  given  on  the  right 
hand  side  of  the  figure.  As  readily  apparent,  even  at  maximum 
current,  the  cell  operating  conditions  are  well  away  from  the 
reactant  starved  regime  where  what  are  normally  termed  con¬ 
centration  losses  become  dominant.  The  streamlines  for  the  flow 
in  the  channels  and  the  electrodes  are  also  plotted  in  Fig.  10  and 
show  the  flow  turning  around  after  hitting  the  electrodes.  The 


Fig.  1 0.  Mass  fraction  profiles  of  the  reactants  in  the  2-D  Maxwell-Stefan  model 
at  0  V,  850  °C. 

streamlines  enter  the  electrolyte  at  the  cathode  and  appear  at  the 
anode  electrolyte  boundary  to  show  how  the  oxygen  is  trans¬ 
ported  across  the  cell  from  the  cathode  through  the  electrolyte 
to  the  anode,  where  it  comes  out  as  H2O. 

5.  Conclusions  and  future  work 

The  primary  ongoing  goal  of  this  work  is  to  build  a  mathe¬ 
matical  model  for  the  single  H2S  SOFC  cell  that  can  then  be 
used  to  help  in  further  development  of  F^S  fuelled  cells  or 
stacks.  Four  models  of  increasing  complexity  have  been  pre¬ 
sented  and  compared  for  a  button  cell  SOFC  fuelled  by  H2S. 
Non-linear  least  squares  were  used  to  estimate  the  unknown 
electrochemical  kinetic  parameters  for  all  the  models.  Exam¬ 
ination  of  these  parameter  estimates,  and  the  fit  between  the 
model  output  and  experimental  data  used,  allows  us  to  identify 
expected  patterns  in  the  parameters  and  compare  the  four  models 
and  their  ability  to  simulate  SOFC  operation.  We  suggest  the  2- 
D  model  with  Maxwell-Stefan  mass  transfer  be  used  in  further 
modelling  studies  because  it  is  the  one  that  accounts  for  the  mass 
transfer  resistance  and  thus  the  activation  resistance  accurately, 
and  is  the  most  comprehensive  of  the  four  models  presented 
here. 

To  comprehensively  validate  these  models,  reliable  estimates 
for  the  exchange  current  densities  from  electrode  characteri¬ 
zation  studies  are  needed.  Detailed  electrochemical  studies  of 
the  electrodes  should  also  throw  some  light  on  the  actual  reac¬ 
tion  mechanisms  and  hopefully  yield  better  rate  expressions.  In 
fact,  given  enough  information  about  the  experimental  setup,  our 
models  can  be  modified  to  simulate  the  electrochemical  charac¬ 
terization  experiments. 

We  believe  that  the  biggest  shortcoming  of  our  model  is  that 
the  reaction  mechanisms  used  to  derive  the  rate  expressions  are 
too  simplistic.  One  of  the  next  steps  in  our  work  will  be  to 
propose  more  complex  reaction  models  that  include  multiple  ele¬ 
mentary  reaction  steps  at  each  electrode  as  well  as  competitive 
adsorption  and  desorption  of  reaction  species.  The  rate  expres- 
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sions  derived  for  the  above  mechanisms  can  then  be  included  in 
our  models  and  should  lead  to  better  matches  to  experimental 
cell  performance  data.  These  rate  equations,  however,  typically 
have  more  kinetic  parameters  than  just  the  exchange  current 
densities.  To  get  reliable  estimates  for  these  kinetic  parameters 
and  to  validate  any  proposed  reaction  mechanisms,  very  careful 
experimental  characterization  of  the  electrodes  using  half  cell 
experiments  [34,35]  is  required. 

Another  possible  approach  to  get  at  electrode  reaction  mecha¬ 
nisms  and  rate  equations  is  through  molecular  modelling.  Meth¬ 
ods  from  theoretical  chemistry  such  as  density  functional  theory, 
which  are  based  on  a  quantum  mechanical  description  of  the 
bonds  between  atoms,  can  be  used  to  find  feasible  electrode 
reaction  mechanisms  [36,37].  These  methods  can  also  give  esti¬ 
mates  for  the  kinetic  parameters  for  the  elementary  reactions  in 
a  given  reaction  mechanism  [38].  Sun  et  al.  [39]  have  used  the 
above  approach  to  investigate  the  hydrodesulfurization  (HDS) 
reaction  on  M0S2,  NiMoS,  and  C0M0S.  In  the  HDS  reaction, 
H2  is  used  to  remove  S  from  the  sulfided  surface,  while  in  a  H2S 
SOFC  anode,  H2S  adsorbs  and  dissociates  on  the  surface.  The 
H  (and  perhaps  the  S)  atoms  near  the  anode-electrolyte  inter¬ 
face  then  presumably  react  with  oxide  ions  in  the  electrolyte  to 
give  reaction  products.  In  future  work,  we  propose  to  use  results 
from  the  approach  outlined  here  to  build  a  better  model  for  the 
electrochemistry  in  the  anode. 

These  updated  models  will  then  be  used  to  do  parametric 
studies  with  variation  in  fuel  composition,  flow-rates  and  tem¬ 
peratures.  We  also  intend  to  include  heat  transfer  in  our  model 
to  examine  the  temperature  variations  in  the  operating  cell  and 
to  check  whether  a  thermally  self  sustaining  SOFC  using  H2S 
is  possible. 
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